#!/usr/bin/python

##
## Copyright (C) 2009 Alberto Ferrarini <alberto.ferrarini@sci.univr.it>
## Copyright (C) 2010 Luciano Xumerle <luciano.xumerle@gmail.com>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.
##

from expression import *

import sys
import os

###########################################################################
## GLOBAL VARS
###########################################################################

alt_splices=dict()
normal_splices=dict()

###########################################################################

def updateCounted( ucSpliceCode ):
   global alt_splices
   global normal_splices
   rr=ucSpliceCode.split("|")
   st=int(rr[3])
   en=int(rr[4])

   if abs(st-en)==1:
       if rr[0] not in normal_splices:
           normal_splices[rr[0]]=dict()
       if (st, en) not in normal_splices[rr[0]]:
           normal_splices[ rr[0] ][ (st, en) ]=0
       normal_splices[ rr[0] ][ (st, en) ]+=1
   else:
       if rr[0] not in alt_splices:
           alt_splices[rr[0]]=dict()
       if (st, en) not in alt_splices[rr[0]]:
           alt_splices[ rr[0] ][ (st, en) ]=0
       alt_splices[ rr[0] ][ (st, en) ]+=1
   return

###########################################################################

def readAlignedFile(vmfFile):
   onlyMulti=False
   mapfields=[]
   ## start reading vmf file
   totalLine=0

   ## VMF file
   if "vmf" in vmfFile:
       for line in open(vmfFile):
           mapfields = splitEland2line(line, onlyMulti)
           if len(mapfields)<1:
               continue
           updateCounted( mapfields[0][6] )
           totalLine+=1

   ## MAQ file
   if "maq" in vmfFile:
       for line in open(vmfFile):
           mapfields = splitMaqLine(line)
           if len(mapfields)<1:
               continue
           updateCounted( mapfields[0][6] )
           totalLine+=1

   ## BOWTIE
   if "vmf" not in vmfFile and "maq" not in vmfFile:
       old=''
       for line in open(vmfFile):
           cur=splitBowtieLine(line)
           if old != cur[0]:
               if (len(mapfields)==1 and not onlyMulti) or (len(mapfields)>1 and onlyMulti):
                   updateCounted( mapfields[0][6] )
                   totalLine+=1
               mapfields=[]
               old=cur[0]
           mapfields.append(cur)
       if (len(mapfields)==1 and not onlyMulti) or (len(mapfields)>1 and onlyMulti):
           updateCounted( mapfields[0][6] )
           totalLine+=1

###########################################################################
# print normal_splices & print alt_splices
###########################################################################

def printCounted( destFile, all_splices, min_reads, norm_factor ):
   f = open(destFile, "w+")
   for id in all_splices:
       output = id + "\t"
       total_reads = 0
       write_it = 0
       for splice in all_splices[id].keys():
           if all_splices[id][splice] >= int(min_reads):
               output = output + str(splice[0]) + ":" + str(splice[1]) + " -> " + str(round(all_splices[id][splice])) + "; "
               write_it = 1
           else:
               output = output + "*" + str(splice[0]) + ":" + str(splice[1]) + " -> " + str(round(all_splices[id][splice])) + "; "
               write_it = 1
       if write_it:
           f.write(output + "\n")
       else:
           if total_reads > 0:
               f.write(output + str(round(total_reads)) + "\n")

###########################################################################

def main():
   """main routine"""
   global alt_splices
   global normal_splices

   # pars
   errorMsg="SYNTAX: " + sys.argv[0] + " <readsCoverage> <factor> <vmfFile>\n"

   if len(sys.argv)<4:
       print >>sys.stderr, errorMsg
       return

   readAlignedFile(sys.argv[3])

   nSplices = sys.argv[3] + "." + sys.argv[1] + ".normsplices"
   aSplices = sys.argv[3] + "." + sys.argv[1] + ".altsplices"

   printCounted( nSplices, normal_splices, sys.argv[1], sys.argv[2] )
   printCounted( aSplices, alt_splices, sys.argv[1], sys.argv[2] )

if __name__=='__main__': main()

